/*	Interactive Gravitational Simulator
	Project: IGS Core
	File: Octree.cpp
	Author: Mike Bantegui <mbante2@gmail.com>, Hofstra University 
	Copyright (C) 2012 - 2013 Mike Bantegui

	This file is part of the IGS software package for simulating
	the N-Body problem in real-time. This was developed as part
	of a departmental honors thesis project at Hofstra University
	in the Spring 2012 semester.

	IGS is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2 of the License, or
	(at your option) any later version.

	IGS is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

#include "Body.h"
#include "Node.h"
#include "Octree.h"
#include <omp.h>

namespace IGS
{

// Constructs an empty octree with the specified opening
Octree::Octree(const double theta) : Tree(theta, Octree::LeafSize)
{
}

// Constructs an empty octree with the specified opening and leaf size
Octree::Octree(const double theta, const unsigned leafSize) : Tree(theta, Octree::BranchSize, leafSize)
{
}

// Destructs the Octree
Octree::~Octree()
{
}

// Returns the name of the Octree
std::string Octree::Name() const
{
	return "octree";
}

#ifdef PoolAllocator
// Calculates the amount of nodes required for a system of N Bodies
unsigned Octree::CalculatePoolSize(const int N) const
{
	return 10 * N / 7;
}
#endif

#ifdef BalancedOctree
Vector CenterOfMass(const std::vector<Body *> &bodies)
{
	double tm = 0;
	Vector mr(0, 0, 0);

	for (size_t i = 0; i < bodies.size(); i++)
	{
		tm += bodies[i]->Mass;
		mr += bodies[i]->Mass * bodies[i]->Position;
	}

	return mr / tm;
}
#endif

// Branches the Octree at the given node, inserting all the bodies as children of the node.
Node *Octree::Branch(std::vector<Body *> bodies, const Vector &min,
	const Vector &max, const unsigned depth)
{
	// Create empty leaf node
	if (bodies.empty())
		return NULL;

	// Create leaf node
	if (bodies.size() <= LeafSize)
		return AllocateNode(bodies);

	const size_t n = bodies.size() / BranchSize;

#ifdef BalancedOctree
	const Vector &center = CenterOfMass(bodies);
#else
	const Vector &center = (min + max) / 2;
#endif

	// Note these are labeled using X-Y-Z notation such as Left-Bottom-Back
	std::vector<Body *> LBB, LBF, LTB, LTF, RBB, RBF, RTB, RTF;
	LBB.reserve(n); LBF.reserve(n); LTB.reserve(n); LTF.reserve(n);
	RBB.reserve(n); RBF.reserve(n); RTB.reserve(n); RTF.reserve(n);

	// Collect the bodies to be inserted in each octant
	for (size_t i = 0; i < bodies.size(); i++)
	{
		Body *b = bodies[i];
		const Vector &r = b->Position;
		if (r.x <= center.x) // Left
		{
			if (r.y <= center.y) // Bottom
			{
				if (r.z <= center.z) LBB.push_back(b); // Back
				else LBF.push_back(b); // Front
			}

			else // Top
			{
				if (r.z <= center.z) LTB.push_back(b); // Back
				else LTF.push_back(b); // Front
			}
		}

		else // Right
		{
			if (r.y <= center.y) // Bottom
			{
				if (r.z <= center.z) RBB.push_back(b); // Back
				else RBF.push_back(b); // Front
			}

			else // Top
			{
				if (r.z <= center.z) RTB.push_back(b); // Back
				else RTF.push_back(b); // Front
			}
		}
	}

	// Free memory being used
	bodies.clear(); bodies.shrink_to_fit();
	std::vector<Node *> branches(Octree::BranchSize, NULL);

	// Recursively branch and insert the bodies into the octants
	#pragma omp parallel sections
	{
		#pragma omp section
		branches[0] = Branch(LBB, Vector(min.x, min.y, min.z), Vector(center.x, center.y, center.z), depth + 1);
		#pragma omp section
		branches[1] = Branch(LBF, Vector(min.x, min.y, center.z), Vector(center.x, center.y, max.z), depth + 1);
		#pragma omp section
		branches[2] = Branch(LTB, Vector(min.x, center.y, min.z), Vector(center.x, max.y, center.z), depth + 1);
		#pragma omp section
		branches[3] = Branch(LTF, Vector(min.x, center.y, center.z), Vector(center.x, max.y, max.z), depth + 1);
		#pragma omp section
		branches[4] = Branch(RBB, Vector(center.x, min.y, min.z), Vector(max.x, center.y, center.z), depth + 1);
		#pragma omp section
		branches[5] = Branch(RBF, Vector(center.x, min.y, center.z), Vector(max.x, center.y, max.z), depth + 1);
		#pragma omp section
		branches[6] = Branch(RTB, Vector(center.x, center.y, min.z), Vector(max.x, max.y, center.z), depth + 1);
		#pragma omp section
		branches[7] = Branch(RTF, Vector(center.x, center.y, center.z), Vector(max.x, max.y, max.z), depth + 1);
	}

	return AllocateNode(branches);
}

}